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Abstract 

^""^ Using simulations of geosynchrotron radiation from extensive air showers, we present a relation between the shape of the geosynchrotron 
i | radiation front and the distance of the observer to the maximum of the air shower. By analyzing the relative arrival times of radio pulses at 
several radio antennas in an air shower array, this relation may be employed to estimate the depth of maximum of an extensive air shower 
if i ts impact position is known, allowing an estimate for the primary particle's species. Vice versa, the relation provides an estimate for the 
. impact position of the shower's core if an external estimate of the depth of maximum is available. In realistic circumstances, the method 
delivers reconstruction uncertainties down to 30 g/cm 2 when the distance to the shower core does not exceed 7 km. The method requires that 
C^he arrival direction is known with high precision. 
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1. Introduction 

o 

One of the most important open questions in astroparticle 
physics is the nature of cosmic-ray particles at the highest ener- 
gies. At energies exceeding 10 15 eV, at present, the only practi- 
cal way to investigate cosmic-ray particles is to register exten- 
ts sive air showers induced by cosmic rays in the atmosphere. In 
t-H such experiments it is only possible to make statements on the 
composition of primary cosmic rays based on statistical eval- 
v~j uations. Abundances of primary particle types of an ensemble 
rN of air showers are frequently derived by looking at the depth 
of the shower maximum, i.e. the depth at which the number of 
particles in a shower reaches its maximum. 

In recent years, there has been a surge of interest in the de- 
tection of extensive air showers by means of the radio emis- 
sion produced by the shower particles [1,2]. This observational 
technique allows one to look all the way up to the shower 
maximum, and it has the advantage over detecting the parti- 
cles themselves at ground level that there is no absorption of 
the signal. Several theories explaining the emission mechanism 
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have been proposed [3, 4, 5]. The former of these explains the 
observed radio emission from the principle of geosynchrotron 
radiation, and using a sophisticated model of geosynchrotron 
emission it was shown that the position of the maximum of 
inclined showers can be derived from the lateral slope of the 
electric field strength at ground level [6]. 

In this work, we use simulations of air showers and their 
geosynchrotron radiation to estimate the value of the depth of 
maximum and the impact position of the shower core. The 
method developed exploits delays in the arrival time of the 
signal at different positions on the ground. 

2. Method 

Detailed distributions of electrons and positrons at differ- 
ent atmospheric depths were obtained from an air shower li- 
brary [7] produced with corsika simulations [8] and the coast 
library [9]. The library contains air showers initiated by pho- 
tons, protons, and iron nuclei of energies in the range 10 16 to 
10 20 5 eV, incident from zenith angles up to 60°. 

A subset of ~ 700 simulations from this library, chosen at 
random, was used to calculate the radio signal emitted by these 
airs showers. The reas code version 2.58 [10, 11] was used to 
obtain the radio pulses associated with each air shower simu- 
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Fig. 1. Layout of the virtual array of radio antennas used in the simulations 
presented in this work. Each marker represents an antenna position. 

lation at an altitude of 100 m above sea level. The expected 
radio signal was calculated for an array of antennas as shown 
in Fig. 1. 

The magnetic field in all simulations, both corsika and reas, 
was taken to match values in northwestern Europe at a field 
strength of 49 /ST and a declination of 68°. The height of the 
detector array was fixed at 100 m above sea level, corresponding 
to an atmospheric depth of X 1024 g/cm 2 . 

3. Parameterization 

For showers hitting the detector at an angle, one has to com- 
pensate for projection effects. Let 0q and <po be the zenith and 
azimuth angle at which the primary enters the atmosphere. For a 
radio antenna a distance d on the ground away from the shower 
core in the direction 6 with respect to the incidence angle <po, 
the impact parameter r is 



cos 2 8 sin 2 



(1) 



The delay r, converted to length units by multiplying with the 
speed of light in vacuum, is defined as the lag of the peak 
strength of the radio signal with respect to the arrival time at 
the shower impact location. It can be written as 



r = t + d cos 5 sin ( 



(2) 



where t(r, 6) is the delay caused by the non-planar shape of the 
shower front expressed in length units. In the analysis in the 
remainder of this work, these geometrical compensations have 
been included. 

In the case of a spherical shower particle front, the expected 
shape of its emitted radio signal is a spherical wavefront as 
well. The delay t can then be written in terms of the distance 
to the center of the sphere R and the distance from the shower 
core r as 

2R' 
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where the approximation holds for r <K R. It was shown previ- 
ously, however, that the assumption of a spherical shower par- 
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Fig. 2. Radio signal delay for a typical vertical 10 185 eV proton shower 
PGnax - 780 g/cm 2 ). Solid curves represent signal delays (converted to length 
units) r at intervals of 5 m (thick lines every 10 m). For reference, perfect 
circles at different distances are also drawn (dotted). 

tide front is unrealistic for large air showers [12]. Therefore, 
the shape of t as a function of r is expected to be different, too. 

The delay of a radio pulse t is defined as the lag between 
a hypothetical plane wave and the actual maximum of the re- 
ceived signal. Fig. 2 shows a contour plot of the distribution 
on the ground of this lag for a typical vertical proton shower at 
E = 10 18 - 5 eV, with X max 780 g/cm 2 . The geomagnetic field 
points north in this figure. Notice the deviation from circularity 
of the front, which is strongest near the shower core in the east 
and west directions. This asymmetry results only from radiation 
processes and is not a consequence of asymmetries in the parti- 
cle front of the shower, because the distributions used to create 
the radio shape are cylindrically symmetric by design [7]. 

Analysis of a set of ~ 700 showers from photons, protons, 
and iron nuclei at various energies and incidence angles as 
described in section 2 reveals that, to first order approximation, 
these delays can be described by the parameterization 



(4) 



where R represents the distance of the impact location to the 
shower maximum, which can be translated unambiguously to 
a value of X max . R] is a scale parameter, the exponent of which 
was chosen to match the dimension of t (distance). Optimisation 
of the parameters reveals that the minimum for the Rq parameter 
is very shallow, and the parameterization can be made to work 
with a wide range of values for it without appreciable change in 
quality of the resulting fit. Therefore, Rq was kept at a constant 
value of 6 km in the final determination of the other parameters 
to speed up the fit process. 

The parameters in the above relation do not depend signif- 
icantly on either primary energy or zenith angle other than 
through the respective influences on the depth of the shower 
maximum. This is not very surprising, because the particle dis- 
tributions responsible for the radiation do not exhibit any de- 
pendence on these parameters either [13, 12]. Though the val- 
ues for Ro, a, and ft depend on the orientation of the shower 
with respect to the magnetic field, this dependence is much 
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Fig. 3. Example of the parameterization presented in (4) and (5) for the 
signal lag for a vertical proton shower at an energy of 10 20 eV and 
Xmax - 895 g/cm 2 . The simulated lag at 6 = 0° and 6 = 90° is indicated by 
crosses and diamonds, and their respective corresponding parameterizations 
are drawn as solid and dashed lines. 

smaller than the average statistical variation between showers. 
Therefore, we will restrict the variations in the parameters to 
a dependence on the angle S only. A fit to the simulated pulse 
lags in the region 40 m < d < 750 m yields the following over- 
all best-fit parameters: 

R\ = 3.87 + 1.56 cos(2(5)+ 0.56 cos 6 (in km), 
a= 1.83 + 0.077 cos(2(5) + 0.018 cos 5, (5) 
P = -0.76 + 0.062 cos(2c5) + 0.028 cos 5. 

The cos(2(5) terms in these equations reflect the asymmetries in 
the east-west versus north-south direction. Note that a < 2 for 
all 5, confirming the non-spherical shape of the wave front. An 
example of the parameterization is shown in Fig. 3, in which 
the simulated lags and their corresponding parameterizations 
are drawn for a vertical proton shower at 10 20 eV and X max ^ 
895 g/cm 2 as a function of distance from the shower impact 
location. Two sets are shown, for 5-0° and 6 = 90°, respec- 
tively. 

The intrinsic accuracy without external error sources of our 
parameterization may be assessed from Fig. 4. This plot shows 
how the distance to the shower maximum R as reconstructed 
from the parameterization in (4) and (5) compares to the actual 
distance as a function of the delay. Note that the figure shows 
reconstructions of single antennas rather than complete show- 
ers: this means that the histogram in this figure is composed of 
80 antennas x 700 showers = 5.6 • 10 4 individual reconstruc- 
tions. It is no surprise that antennas with longer delays of t > 
10 m produce more accurate reconstructions, since the relative 
error is smaller there. Even at arrival lags of less than 1 m, 
however, the standard deviation is less than 10 % of the actual 
value. 

In a typical array of radio antennas, one can determine the 
delays r very accurately: using modern equipment, resolutions 
down to a few ns can be achieved. We can use the delay val- 
ues to employ the parameterization in (4) in two ways: if the 
position of the shower core is known accurately by scintilla- 
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Fig. 4. Relative intrinsic error in the reconstruction of R as a function of 
the delay t. Darker areas mark higher numbers of reconstructions. The total 
amount of colouring is constant for every slice in r; the intensity is in arbitrary 
units. 

tor measurements, we can use it to estimate the distance to the 
shower maximum. If, on the other hand, an estimate for the 
depth of maximum is available, the position of the shower core 
can be reconstructed. We explore these possibilities in detail in 
the following two sections. 

4. Determining depth of shower maximum 



By rearranging (4), we may write 



R = R 



1 \f) 



■Ro 



(6) 



to reconstruct the distance to the shower maximum. Using this 
parameterization, the reconstructed distance to the shower max- 
imum is plotted versus the simulated value in the left panel of 
Fig. 5. Each dot in this plot represents the reconstructed value 
of R for one shower event, obtained by taking a weighted aver- 
age of the reconstructions from the delays in individual anten- 
nas. If the antennas are placed on a regular grid, a weight oc r 2 
seems justified to match each time delay to its expected rela- 
tive error, since a 2. Our simulated array is denser near the 
shower core, which was compensated for by multiplying by an 
extra factor of r, arriving at a total weight for each antenna oc r 5 . 

Around each mark in Fig. 5 a circle is drawn, the radius of 
which is the distance corresponding to an atmospheric depth 
of 20 g/cm 2 at the position of the simulated air shower maxi- 
mum. This value represents the average error for reconstructed 
X max values with the Pierre Auger Observatory using air fluo- 
rescence techniques [14]. The algorithm correctly reconstructs 
the distance to the shower maximum as simulated, with a stan- 
dard deviation of 216 m. Note that for negative distances, the 
shower maximum lies below the observation level. By design 
of he algorithm, correct reconstruction of these negative dis- 
tances is possible, but only if the downward distance is smaller 
than Rq. 

When the uncertainties in Fig. 5 are converted to atmospheric 
depths, we find that the standard deviation of the values for 
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Fig. 5. Scatter plot for ~ 700 showers of various species and energies E > 10 17 eV of simulated values for R versus the values as reconstructed by the method 
outlined in the text. Circles around each reconstruction represent error margins of 20 g/cm 2 . The left panel shows the theoretical limit in reconstruction 
accuracy, while in the right plot realistic Gaussian errors were introduced around the observables in (6). Note that the distance to the shower maximum extends 
below zero: these are showers that reach their maximum below the level of the observing radio array. 



AX max is between 15 and 20 g/cm 2 over the full energy range 
of 10 16 -10 2l) eV. However, we have so far considered perfect 
circumstances, assuming exact knowledge of the impact angle 
and position of the shower axis as well as the delay of the radio 
pulses. A more realistic picture emerges by introducing some 
error sources in the reconstruction. For a dense array of radio 
antennas, such as the lopes [1] or lofar [15] telescopes, the ac- 
curacy in the arrival direction is of the order of 1.0° [16]. A fea- 
sible time resolution for determining the maximum pulse height 
is about 10 ns (3 m). Because errors in the antenna positions 
can be reduced to less than 10 cm by extended gps measure- 
ments, they do not contribute significantly to this uncertainty. 
The accuracy in determining the position of the shower core 
has not been investigated thoroughly yet using radio detection. 
Therefore, we adopt a typical value from the analysis of the 
kascade experiment data of 1 m [17, 18]. All of the above er- 
rors are assumed to follow Gaussian distributions. Additionally, 
we ensure that the signal is sufficiently strong by demanding a 
field strength over 180 //V/m, which corresponds to a power 
signal-to-noise ratio of around 3 in a rural area [6]. 

The right panel of Fig. 5 shows the situation when these er- 
ror estimates are included. The correlation is reduced signifi- 
cantly, which is mainly the result of the uncertainty in the ar- 
rival direction of the shower. For very inclined showers in par- 
ticular this can change the expected delay times dramatically. 
When the accuracy of the shower impact location is reduced, 
this mostly affects showers for which the maximum lies at a 
large distance from the observer. When the error is increased 
to 5 m, for example, hardly any predictions can be made for 
distances > 10 km. 
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Fig. 6. Distribution of residuals for the reconstruction of the depth of maxi- 
mum for various primary energies. Plots are shown for urban, rural, and ideal 
noise level scenarios. 

The distribution of residuals AX nYdx (i.e. the reconstructed 
minus the simulated value of the depth of maximum) is shown 
in Fig. 6 for primary energies between 10 17 and 10 20 eV. In this 
plot, a homogeneous detector sensitivity up to zenith angles 
8 < 60° is assumed. Three background noise scenarios are 
shown: one for an ideal noise level (requiring a field strength 
|£1 > 65 /iV/m for successful determination of t), one for a 
rural environment (|£] > 180 fiV/m), and one corresponding 
to an urban area QE\ > 450 yuV/m) [6]. 

From this figure, we observe that the reconstruction accu- 
racy for X max decreases rapidly at low energies. This is because 
low-energy showers do not occur very deep in the atmosphere 
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Fig. 8. Density plot for ~ 700 showers of various species and energies E > 10 17 eV of simulated values for the impact location of the shower as reconstructed 
by the method outlined in the text. The actual position of the core is marked with a cross. Also shown is the arrival direction for slanted air showers (solid 
line). The left panel shows the theoretical limit in reconstruction accuracy. On the right realistic observational errors were introduced in (7). The colour 
intensity scales linearly with the number of reconstructions at that point. 
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constant, and the background noise was 65 //V/m. The values 
at 1 m correspond to the distribution widths in the rightmost 
panel in Fig. 6. From Fig. 7 it is clear that the reconstruction 
technique employed in this section is better suited for dense 
arrays, where more accurate impact locations are available. 

If the maximum available distance to the shower core is very 
small, as would be the case for an array such as lopes, the 
fraction of good reconstructions is reduced dramatically. This 
makes sense, as the shower front shape can no longer be probed 
accurately. In particular, if the radius of the array decreases 
to less than ~ 500 m, the amount of useful reconstructions is 
negligible. 

5. Determining shower core position 



Fig. 7. Dependence of cr(AX mllx ) on the uncertainty introduced in the air 
shower's impact location for showers of different energies. 



on average, raising the distance to the shower maximum, espe- 
cially in slanted showers. This results in a radiation front with 
less curvature, necessitating delay measurements further away 
from the impact location to obtain the same level of reconstruc- 
tion accuracy. The produced field strength, however, is propor- 
tional to the primary energy, decreasing the patch size that is 
sufficiently illuminated. The combined effect is that it is hard 
to make correct estimations for the depth of maximum of low 
energy showers, unless an array at high altitude is employed. 

Additionally, the behaviour of the reconstruction accuracy 
curve at 10 18 eV in the three scenarios highlights the impor- 
tance of low background interference levels: the width of the 
distribution decreases dramatically at this energy. It is also ob- 
served that the distribution width does not vary much for ener- 
gies of 10 19 f and 10 20 eV. 

Fig. 7 shows the resulting standard deviation in the values for 
AX max when the uncertainty in the impact location of the shower 
is varied. Gaussian error values on other parameters were held 



If an estimate for X max (and therefore for R) is available, we 
can employ (4) in an alternative way to estimate values for the 
impact parameter r, by writing 



R 



l+l/a/3-l/a 
1 



Ala 



(7) 



(R + RoyW 

In an actual experimental setting, the dependencies of a, 
ft, and Ri on 5 need to be taken into account, for example 
through an iterative fitting procedure for r and 6. For the sake 
of simplicity, we will only reconstruct the distance to each 
antenna here, and we will assume the general direction of the 
core impact position to be known. This decision is motivated 
by the fact that the effect on the value of r caused by variations 
in 5 is generally small. 

In the theoretical limit, the distribution of reconstructed 
shower core positions using this method is shown in the left 
panel of Fig. 8. The colouring in this plot shows the amount of 
reconstructions at a certain position relative to the actual core 
impact location. The true position is at the origin, indicated 
by a cross. The arrival direction of inclined showers is always 
from the left, as indicated by the arrow. Note that the elongated 
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Fig. 9. Effects on the maximum field strength to the east of a vertical 10 18 eV 
shower arising from applying a rectangular filter to the raw radio pulse. From 
top to bottom, the unfiltered pulse is shown, and the same pulse with filters 
of 10, 20, 30, and 40-100 MHz applied. 

structure of the reconstruction distribution is not a projection 
effect from inclined showers: we have already compensated 
for this by the transformation to the shower plane through (1). 
Instead, the feature is a systematic error intrinsic to the re- 
construction algorithm. For a shower incident from the south, 
for example, the parameterized form is not symmetric in the 
north-south direction, but it is in the east-west direction. This 
effect is also responsible for the slight offset of nearly -2 m in 
the x direction. 

Theoretically, the systematic offset could be reduced and 
possibly even removed entirely by refining the parameterization 
in (4) and (5). There is little gain in this exercise, however, when 
a more realistic reconstruction estimate is made. This is clarified 
in the right panel of Fig. 8, where again some error sources were 
introduced. The error in the arrival direction is again 1.0°, and 
a Gaussian uncertainty of 20 g/cm 2 in the value of the shower 
maximum is assumed, corresponding to a typical error in R 
of 200-250 m. Clearly, the offset mentioned earlier is entirely 
swamped by the deviations induced by the uncertainties. The 
substantial difference in reconstruction accuracy between the 
x and y direction results directly from the uncertainty imposed 
on 9q : even a small deviation of the zenith angle will make a 
noticeable difference in the obtained value for t from (2). 

Similar to the determination of X max , the average error in- 
creases drastically when the radius of the array is smaller than 
500 m. The error does not increase significantly, however, when 
the minimum distance is set to 300 m. This is slightly counter- 
intuitive, but it is again related to the accurate probing of the 
shower front shape. Of course, the requirement remains that the 
arrival delay at the impact location is known to 10 ns or so. 

6. Discussion 

The analysis in this work on the relative delays of geosyn- 
chrotron emission from extensive air showers was performed 
on the raw, unfiltered pulse shape. In real experiments, how- 
ever, the antennas used are bandwidth-limited, which will be 



reflected in the measured field strength of the pulse. This is 
illustrated in Fig. 9, which shows the maximum value of the 
measured field strength when different rectangular filters are 
applied to the raw pulses. When frequencies below 40 MHz 
are clipped, the field strength is around 10% of the unfiltered 
value over the entire distance range. 

Another effect that has not been investigated is that of the 
observer's altitude: in our simulations, this height was fixed at 
100 m above sea level. We do not anticipate a significant change 
of the parameterization or its parameters, however. This can be 
inferred from the fact that the description is valid independent 
of zenith angle. Changing this angle is comparable to varying 
the observer's altitude. 

Though a deviation from a planar wave is indeed observed in 
lopes measurements [1], at only 200 m the array is too small to 
benefit from the theoretical knowledge of the shape of the radio 
pulse front. There are currently two other experiments under 
construction, however, that could make use of the technique 
outlined in this work. One of these is the initiative in which 
radio antennas inside the Pierre Auger Observatory [19] will 
be erected [20]. Such an array could use the method in Sect. 5 
to increase the accuracy of the estimated core impact position, 
since its reconstruction error for the surface detectors is in 
excess of 100 m. A precise estimate for X max would have to be 
provided by the fluorescence detectors. The planned spacing of 
radio antennas is 150-375 m, which would allow an accuracy 
in the reconstruction of around 30 m if the core lies within the 
radio array. Using the Auger array for the metod outlined in 
Sect. 4 would probably not be possible, as the uncertainty in 
the reconstructed core position of around 150 m would wash 
out any sensitivity of the algorithm to the shower maximum. 

Another possible experiment is the lofar telescope [15], 
which consists of a dense core of approximately 2 km in di- 
ameter, with groups of 48 radio antennas every few hundred 
meters. Its size and spacing make this setup ideally suited to de- 
termine X m . dx using the method outlined in Sect. 4. The shower 
core position, which would have to be known to apply the 
method, could be obtained in several ways. First of all, there is 
a small scintillator array coincident with the lofar core, allow- 
ing an independent measurement of this quantity. Alternatively, 
pulse shape and lateral slope of the radio signal could be used 
to get an estimate for the core position [21]. It is assumed that 
reconstruction with a dense radio array such as lofar, which 
places antennas at distances of the order of 10 m, will be on a 
par with the precision level of scintillator arrays. 

Pulse shape and lateral slope also contain additional informa- 
tion about the value of X max , with precisions of up to 16 g/cm 2 
[6]. Ideally, one would combine the two methods in a single fit 
to obtain the best possible reconstruction accuracy. 

7. Conclusion 

Through detailed simulations of air showers and their 
geosynchrotron radio emission, we have derived an empirical 
relation between the relative delay of the radio pulse emit- 
ted by the air shower front and the atmospheric depth of the 
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shower maximum. By analysis of the radio pulse arrival delays 
in radio antennas in an array of low-frequency radio antennas, 
this relation can be used to estimate the depth-of-maximum if 
the impact position is known or vice versa. 

We have confirmed that both methods work in principle, with 
no information other than radio signal delays used in the re- 
construction. When the algorithm is tested under realistic con- 
ditions, however, the accuracy of the method is reduced. In the 
case of determining the shower maximum, reconstruction down 
to a useful confidence level is possible only for shower max- 
ima up to ~ 7 km away, and only if the shower core impact 
position is known down to a few meters. When the parameteri- 
zation is used to derive this position, the critical quantity is the 
accuracy in the zenith angle of the shower, which needs to be 
significantly less than a degree to reconstruct the shower impact 
location to an accuracy of 10 m at high inclinations up to 60°. 
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